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Abstract.  A  discontinuous  Galerkin  (DG)  finite  element  formulation  is  proposed  for  the  solu¬ 
tion  of  the  compressible  Navier-Stokes  equations  for  a  vertically  stratified  fluid,  which  are  of  interest 
in  mesoscale  nonhydrostatic  atmospheric  modeling.  The  resulting  scheme  naturally  ensures  con¬ 
servation  of  mass,  momentum,  and  energy.  A  semi-implicit  time-integration  approach  is  adopted 
to  improve  the  efficiency  of  the  scheme  with  respect  to  the  explicit  Runge-Kutta  time  integration 
strategies  usually  employed  in  the  context  of  DG  formulations.  A  method  is  also  presented  to  refor¬ 
mulate  the  resulting  linear  system  as  a  pseudo-Helmholtz  problem.  In  doing  this,  we  obtain  a  DG 
discretization  closely  related  to  those  proposed  for  the  solution  of  elliptic  problems,  and  we  show 
how  to  take  advantage  of  the  numerical  integration  rules  (required  in  all  DG  methods  for  the  area 
and  flux  integrals)  to  increase  the  efficiency  of  the  solution  algorithm.  The  resulting  numerical  for¬ 
mulation  is  then  validated  on  a  collection  of  classical  two-dimensional  test  cases,  including  density 
driven  flows  and  mountain  wave  simulations.  The  performance  analysis  shows  that  the  semi-implicit 
method  is,  indeed,  superior  to  explicit  methods  and  that  the  pseudo-Helmholtz  formulation  yields 
further  efficiency  improvements. 
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1.  Introduction.  In  recent  years,  great  attention  has  been  devoted  to  the  dis¬ 
continuous  Galerkin  (DG)  finite  element  method  in  the  context  of  geophysical  fluid 
dynamics  applications.  This  is  motivated  by  the  fact  that  the  DG  framework  simul¬ 
taneously  provides  a  high-order  discretization,  great  flexibility  in  the  choice  of  the 
computational  grid,  discrete  balance  relations,  robustness  with  respect  to  unphysical 
oscillations,  and  compact  computational  stencils  which  are  a  key  element  in  order  to 
exploit  distributed-memory  parallel  computers  with  up  to  tens  of  thousands  of  proces¬ 
sors.  Without  attempting  to  provide  a  complete  review  of  the  literature,  we  mention 
here  [47,  2,  30,  39,  34,  32],  where  DG  shallow  water  models  are  presented.  The  ap¬ 
plication  of  the  DG  method  to  compressible,  nonhydrostatic  atmospheric  flows,  using 
the  Navier-Stokes  equations  or,  when  the  flow  is  assumed  to  be  inviscid,  the  Euler 
equations,  is  then  considered  in  [31],  where  it  is  shown  that  the  method  represents 
a  good  candidate  for  the  development  of  numerical  climate  and  weather  models.  In 
the  present  paper,  we  continue  the  study  initiated  in  [31]  by  focusing  on  the  aspect 
of  the  time  discretization  which  is,  in  fact,  the  most  penalizing  drawback  of  the  DG 
method  due  to  its  high  computational  cost.  This  latter  cost  stems  from  the  following 
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reasons:  the  need  of  a  complete  independent  set  of  points  for  each  element  to  span 
the  local  polynomial  space;  a  stricter  time-step  than  in  continuous  methods  due  to 
the  upwinding  mechanism  in  the  numerical  flux;  and  the  fact  that  all  previous  ap¬ 
plications  have  been  confined  to  explicit  time-integration.  Our  goal  is  to  remedy  the 
last  two  points.  More  precisely,  the  main  purpose  of  this  article  is  to  show  that  it 
is  possible  to  use  a  semi-implicit  (SI)  time-stepping  method  in  combination  with  a 
DG  spatial  discretization,  thus,  improving  the  computational  efficiency  with  respect 
to  the  explicit  time  discretization  approach,  without  affecting  the  positive  features  of 
the  method.  The  idea  of  numerically  integrating  in  time  a  differential  equation  by 
using  an  explicit  treatment  for  the  nonstiff  component  and  an  implicit  treatment  for 
the  stiff  component,  thus,  obtaining  a  method  which  can  be  generally  classified  as 
“semi-implicit”,  is  very  general  and  has  been  exploited  in  different  contexts,  where 
correspondingly  different  criteria  have  been  used  to  identify  the  stiff  terms.  In  the 
context  of  geophysical  fluid  dynamics,  the  SI  approach  was  first  introduced  in  [38] 
for  the  solution  of  the  primitive  hydrostatic  equations,  where  the  stiff  component  was 
identified  with  the  terms  responsible  for  the  gravity  waves.  Modifications  were  pro¬ 
posed  in  [10]  for  the  hydrostatic  case  and  the  technique  was  then  extended  to  the 
complete  nonhydrostatic  Euler  equations  in  [53,  18,  52],  with  implicit  treatment  of 
both  acoustic  and  gravity  waves.  Since  then,  the  method  has  been  widely  employed 
in  climate  and  weather  prediction  models,  and  we  refer  to  [43,  7,  33]  for  further  de¬ 
tails.  To  our  knowledge,  this  type  of  SI  time  discretization  has  never  been  used  in 
combination  with  a  DG  spatial  discretization,  and  the  SI-DG  formulation  described 
in  the  present  paper  represents  a  novel  contribution.  This  formulation,  moreover, 
besides  being  suited  for  atmospheric  flow  problems,  may  be  useful  also  in  the  context 
of  smaller  scale,  low  Mach  number  fluid  dynamics  simulations. 

We  notice  here  that  a  SI  time  discretization  combined  with  a  DG  space  discretiza¬ 
tion  has  been  considered  in  [56,  55]  for  scalar  equations  with  higher  order  derivatives, 
where  the  stiffness  is  due  to  the  presence  of  the  high  order  terms,  and  in  [23,  26,  27] 
for  the  Navier- Stokes  equations.  Concerning  this  latter  formulation,  we  observe  that 
it  is  different  from  the  one  we  propose  in  the  present  article.  The  method  originally 
developed  by  Dolejsi  and  Feistauer  in  [23]  relies  on  a  careful  handling  of  the  nonlinear 
problem  arising  from  a  fully  implicit  time-stepping  strategy,  and  a  key  ingredient  of 
the  method  is  the  choice  of  a  particular  numerical  flux.  This  approach,  thus,  allows 
the  use  of  large  time-steps  regardless  of  the  Mach  number.  On  the  contrary,  in  our 
approach  only  the  terms  responsible  for  acoustic  and  gravity  waves  are  treated  implic¬ 
itly,  and  there  is  much  freedom  in  the  choice  of  the  numerical  flux.  As  a  result,  the 
scheme  we  propose  is  less  restrictive  in  regards  to  the  numerical  flux  and  is  suitable 
for  the  treatment  of  low  Mach  number  flows,  where  it  results  in  a  simpler  implicit 
problem  than  [23]  since  the  advective  contributions  are  absent. 

A  final  comment  is  in  order  concerning  the  solution  algorithm  for  the  linear  sys¬ 
tem  associated  with  the  SI  time  discretization.  In  the  present  paper,  we  show  that  it 
is  possible  to  reformulate  such  a  system  as  a  pseudo-Helmholtz  problem  for  the  sole 
pressure  variable,  the  discretization  of  which  is  closely  related  to  the  local  discontinu¬ 
ous  Galerkin  (LDG)  methods  proposed  in  [13,  3]  for  elliptic  problems.  By  doing  this, 
the  number  of  unknowns  is  reduced  by  a  factor  four  (and  the  resulting  matrix  prob¬ 
lem  by  a  factor  of  16),  since  the  sole  pressure  variable  appears  in  the  linear  system 
rather  than  the  complete  set  of  degrees  of  freedom  composed  by  density,  horizontal 
and  vertical  velocities,  and  total  energy.  Due  to  this  reduction  in  the  size  of  the  im¬ 
plicit  system,  the  reformulation  as  a  pseudo-Helmholtz  problem  yields  a  significant 
efficiency  improvement  compared  to  the  solution  of  the  original  system. 
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An  outline  of  the  paper  is  as  follows.  In  section  2  the  governing  equations  are 
introduced.  Section  3  deals  with  the  SI  time  discretization,  while  the  DG  spatial  dis¬ 
cretization  is  presented  in  section  4.  The  space-time  fully  discretized  problem  is  then 
summarized  in  section  5.  In  particular,  a  method  to  reformulate  the  resulting  linear 
system  as  a  pseudo-Helmholtz  problem  is  described  in  section  5.2.  The  numerical 
validation  of  the  proposed  formulation  is  presented  in  section  6.  Finally,  conclusions 
and  future  developments  are  discussed  in  section  7. 

2.  Governing  equations.  In  this  section,  we  introduce  the  continuous  equa¬ 
tions  representing  the  mathematical  model  for  the  atmospheric  flow  problem  con¬ 
sidered  in  the  present  paper.  Various  alternative  equation  sets  have  been  proposed 
in  the  literature  to  describe  the  nonhydrostatic  flow  of  a  dry,  stratified  atmosphere. 
In  [31]  three  different  equation  sets  are  examined  by  comparing  the  results  of  five 
spectral  element  and  discontinuous  Galerkin  codes,  and  it  is  found  that  the  conserva¬ 
tive  Navier  Stokes  equations  using  density,  momentum,  and  total  energy  as  prognostic 
variables  represent  one  of  the  most  effective  choices  as  far  as  accuracy  is  concerned.  In 
addition,  the  fact  that  such  an  equation  set  is  in  conservation  form,  even  when  taking 
into  account  viscous  stresses,  makes  it  a  suitable  starting  point  for  the  construction  of 
a  numerical  scheme  endowed  with  discrete  conservation  properties.  For  these  reasons, 
this  equation  set  is  considered  in  the  present  work.  Restricting  ourselves  for  simplic¬ 
ity  to  a  two-dimensional  case  in  the  vertical  (x,  z )  plane,  and  neglecting  the  Coriolis 
terms,  we  obtain  the  following  system  (see  [5]): 

(2.1)  ^+V-^e(q)-V-^(q,Vq)  =  G(q), 

where  each  term  is  defined  as  follows:  q  =  (p,  \T ,E)T  are  the  conserved  quantities, 
p  is  the  density,  and,  letting  v  denote  the  velocity  field,  cv  and  g  denote  the  specific 
heat  for  constant  volume  and  the  gravitational  constant,  respectively,  T  denote  the 
temperature  and  e  =  cvT  +  •  v  +  gz  represent  the  total  energy,  V  =  /tv  is  the 

momentum  and  E  =  pe  is  the  energy  density.  The  inviscid  and  viscous  fluxes  and  the 
source  term  in  (2.1)  are  given  by 


r  v  1 

o 

'  o  ] 

•fe(q)  =  < 

±V®V  +pl 

V  i*V  J 

►,  ^(q,Vq)  =  < 

^(q,vq) 

[  F|(q,Vq)  J 

>,  G(q)  =  < 

PS 

{  0  J 

where  p  is  the  pressure,  g  =  —g k,  with  k  =  (0, 1  )T ,  H  =  E  +  p  is  the  enthalpy,  and 
the  viscous  fluxes  are  defined  as 

Fr  =  M  [Vv  +  VvT  +  AV  •  vl] ,  +  v  • 


p  and  A  denoting  the  two  viscosity  coefficients,  cp  denoting  the  specific  heat  for 
constant  pressure,  and  Pr  being  the  Prandtl  number.  Using  the  Stokes  hypothesis, 
we  can  set  A  =  —  |.  Notice  that,  although  for  practical  applications  a  turbulence 
closure  relation  has  to  be  considered  to  define  the  viscosity  coefficients,  in  the  present 
work  we  assume  that  these  latter  are  known  constants.  Closure  of  system  (2.1)  is 
obtained  through  the  equation  of  state,  which,  in  terms  of  the  solution  variables,  is 
written  as 

(2.2)  p=R(E-l-—-pgzY 

cv  \  2  p  J 

with  R  =  cp  —  cv.  Equation  (2.1)  reduces  to  the  Euler  equations  when  p  =  0. 
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When  solving  system  (2.1)  for  atmospheric  flows,  it  can  be  expected  that  the 
flow  is  nearly  hydrostatic,  i.e.,  in  the  vertical  momentum  equation  the  two  terms 
(from  the  inviscicl  flux)  and  pg  (from  the  source  term)  are  much  larger  than  the 
remaining  ones.  This  can  cause  instabilities  in  the  numerical  approximation  of  the 
problem,  due  to  cancellation  of  significant  digits.  To  avoid  this  effect,  problem  (2.1) 
is  usually  reformulated  in  terms  of  deviations  from  a  constant-in-time  reference  state 
(see  [25]  and,  for  an  alternative,  more  sophisticated  approach,  [8]).  Thus,  we  introduce 
q  =  (p,  VT,  E)t  and  p  =  p(p ,  V,  E)  such  that  V  =  0,  J|  =  0,  ff=0  (also  implying 
f  -  0),  and 

p.3)  I  -  -T»- 

Upon  defining  q'  =  q  —  q  and  p'  =  p  —  p,  we  obtain 
8q' 


(2.4) 

where 


dt 


+  V-^c,(q)-V-^’(q,Vq)  =  G(q') 


^e'(q)  = 


iV( 

p 


V 

)V  +p'l 
i  HV 

p 


Equation  (2.2)  allows  for  an  expression  of  p'  which  is  independent  from  p,  E 


R 


1  U2 


P'  =  ~  E'-- - p'gz 


so  that  no  cancellation  problems  occur  in  evaluating  the  pressure  perturbation.  Prob¬ 
lem  (2.4)  has  the  advantage  that  all  the  terms  in  the  vertical  momentum  equations 
are  of  the  same  order  of  magnitude.  For  this  reason,  in  the  following  we  will  always 
use  (2.4)  instead  of  (2.1),  and  we  will  drop  the  primes  for  simplicity. 


2.1.  Treatment  of  the  open  boundary  conditions.  In  practical  applications, 
it  is  usually  necessary  to  truncate  the  computational  domain  with  artificial  boundaries, 
not  corresponding  to  any  physical  entity.  Ideally,  an  “open  boundary”  condition  is 
desired  on  these  boundaries,  avoiding  any  reflection  of  outgoing  signals.  A  simple  and 
robust  solution  is  represented  by  an  absorbing  layer,  also  known  as  a  sponge  layer, 
as  discussed  in  [36,  24].  In  these  references,  a  method  is  described  to  optimize  the 
structure  of  the  absorbing  layer  in  the  case  where  the  primitive  variable  formulation 
is  considered,  i.e.,  when  the  prognostic  variables  are  Exner  pressure,  velocity,  and 
potential  temperature.  To  apply  these  results  to  the  conservative  formulation  (2.4),  we 
can  in  principle  perform  a  change  of  variables  from  conservative  to  primitive  variables, 
apply  the  damping  coefficients,  and  then  transform  back  to  conservative  variables. 
In  practice,  this  can  be  done  as  follows:  denoting  by  qt  =  {pb,Vf ,  Eb)T  a  known 
boundary  datum,  we  modify  system  (2.4)  to  obtain 

(2.5)  ^=^JVS(q)-TCT(q-q6), 

with  yNS  (q)  =  -V  •  ^(q)  +  V  ■  fv(q,  Vq)  +  G(q)  and 

(2.6)  Tcv  =  M~xTpvM. 

In  (2.6),  M  is  a  4  x  4  matrix  representing  the  linearized  transformation  between 
conservative  variables  and  primitive  variables,  the  linearization  being  performed  in  a 
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neighborhood  of  q while  Tpv  is  a  diagonal  matrix  whose  four  entries  are  computed 
as  in  [36,  24].  Once  the  absorbing  layer  has  been  introduced  within  the  compu¬ 
tational  domain,  the  particular  boundary  condition  prescribed  on  the  boundary  of 
the  computational  domain  itself  has  in  practice  no  effect  on  the  computed  solution, 
so  that,  because  of  the  ease  of  implementation,  we  impose  the  Dirichlet  condition 
q  =  qh.  We  notice  that,  although  introducing  an  absorbing  layer  is  a  common  solu¬ 
tion  to  handle  nonreflecting  boundary  conditions,  more  sophisticated  alternatives  are 
possible.  In  particular,  we  are  currently  exploring  Higdon-type  high-order  boundary 
conditions  [20,  22,  21]  which  we  reserve  for  future  work  since  this  topic  is  beyond  the 
scope  of  the  present  paper. 

2.2.  Identification  of  the  fast  waves  in  the  model.  In  order  to  devise  a 
semi-implicit  time-integration  scheme  for  (2.4),  it  is  first  necessary  to  identify  the 
terms  responsible  for  the  fastest  waves  in  the  model,  i.e.,  acoustic  and  gravity  waves. 
In  fact,  in  the  semi-implicit  time-integration  procedure,  an  implicit  treatment  will 
be  selectively  applied  to  such  terms,  while  an  explicit  approach  will  be  adopted  for 
the  remaining  ones.  Following  [37],  we  identify  these  terms  as  the  divergence  of  the 
mass  flux  in  the  continuity  equation,  the  pressure  gradient  and  the  buoyancy  term 
in  the  momentum  equation  and  the  divergence  of  the  enthalpy  flux  in  the  energy 
equation  (see  also  the  simplified  stability  analysis  in  [42]).  To  linearize  these  terms, 
we  introduce  the  linear  operator 

(2.7)  ^NS(q)  =  -V-^(q)  +G(q) 

with 

r  v 

Fse(d)  =  <  £( E-pgz)l 

[  hV 

and  h  =  =  (E  +p).  The  linearized  form  (2.7)  represents  the  basis  for  the  definition 
of  the  linear  problem  which  will  be  discussed  in  section  4.4. 

3.  Temporal  discretization.  In  this  section  we  discuss  the  semi-implicit  meth¬ 
od  adopted  for  the  time  discretization  of  system  (2.5).  We  notice  that,  although  the 
SI  approach  is  usually  considered  in  combination  with  either  the  Crank-Nicolson  or 
the  leapfrog  time  integration  schemes,  following  the  abstract  formulation  of  [33]  it  is 
possible  to  include  other  time-integration  schemes,  such  as  the  backward  difference 
scheme  presented  in  [35]  and  summarized  in  the  sequel.  We  also  notice  that  all  these 
schemes  can  be  modified  to  account  for  variable  time-steps;  however,  we  consider 
here  for  simplicity,  the  constant  time-step  case  and  leave  the  adaptive  case,  with  the 
associated,  nontrivial,  issue  of  defining  a  criterion  for  varying  the  time  increment,  for 
a  future  work.  We,  thus,  consider  the  abstract  problem 

(S-1)  §  =  >(q), 

to  be  solved  in  ( 0,Tf  i„ ]  with  a  suitable  initial  condition,  define  an  affine  operator  srf 
such  that  srfq  «  J^(q)  and  set 

(3.2)  ^  =  {^(q)  -^q}  +^q. 

The  main  idea  is  now  to  treat  the  term  in  braces  explicitly,  while  the  remaining  term 
will  be  treated  implicitly.  Notice  that  problem  (2.5)  can  be  recast  in  the  form  (3.1)  by 
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Table  3.1 

Thed-method  (TM),  leapfrog  (LF2),  and  backward  difference  of  order  2  (BDF2)  and  3  (BDF3) 
time-integration  schemes  with  their  associated  coefficients  in  the  context  o/(3.4). 


Method 

Q!0 

ai 

Ot-2 

7 

do 

01 

02 

(7-1 

^0 

<7 1 

C 72 

TM 

1 

0 

0 

1 

1 

0 

0 

e 

l-6» 

0 

0 

LF2 

0 

1 

0 

2 

1 

0 

0 

e 

0 

1  —  9 

0 

BDF2 

4 

3 

1 

3 

0 

2 

3 

2 

-1 

0 

1 

0 

0 

0 

BDF3 

18 

11 

9 

11 

2 

11 

6 

11 

3 

-3 

1 

1 

0 

0 

0 

setting  f/  (q)  =  SdNS  (q  )—Tcv  (q  —  q&).  To  allow  for  more  flexibility  in  the  treatment 
of  the  damping  term,  it  is  also  convenient  to  further  split  the  affine  operator  srf  as 
+  &/0,  and  set  s/q  =  d£q  +  f  and  £/°q  =  A?°q  +  f°,  where  f  and  f°  are 
assumed  to  be  constant  in  time.  Equation  (3.2)  can  now  be  written  as 

(3.3)  ^  =  {^(q)  -J?q}+J?q  +  J?°q  +  f°, 

at 

where  S*(q)  =  J^{q)  —  &/°q.  Clearly,  for  a  specific  operator  5?  in  (3.1),  a  proper 
definition  of  the  linear  and  nonlinear  operators  in  (3.3)  is  a  critical  step  in  the  design 
of  the  numerical  scheme.  For  the  case  of  the  Navier-Stokes  system,  this  will  be 
discussed  in  detail  in  section  5.  For  the  present  section,  however,  it  is  sufficient  to 
consider  the  abstract  form  (3.3).  Choosing  now  a  time-step  At,  letting  tn  =  nAt, 
with  n  =  0, . . .  ,Tf  in/At ,  and  denoting  by  q"  the  approximate  solution  at  time  level 
tn,  the  discretization  of  (3.3)  is  constructed  as  follows: 


(3.4) 


dq 

dt 


1 

7  A  t 


2 


2  amq‘ 

m= 0 


2 

{y(q)  -X q}  *  2  ft"  (^(q”“m)  - 

m= 0 


2 


m=  —  1 

^°q  «  ffff°qn+1 


for  suitable  coefficients  7,  am,  /3m,  and  am.  Notice  that,  for  consistency,  it  is  required 
that  am  =  Ylm  pm  =  Yjm  <jm  =  1.  Typically,  the  operator  2z?  is  chosen  in  such  a 
way  that,  for  a  particular  range  of  q,  the  term  Sf — 2zf  vanishes,  and  time  integration  is 
performed  with  the  implicit  scheme  (3.4)  1^4.  In  Table  3.1  it  is  shown  how  to  recover 
some  classical  time-marching  schemes  by  properly  choosing  the  coefficients  in  (3.4), 
namely,  the  0-method,  the  leap-frog  method,  the  backward  difference  scheme  (BDF) 
proposed  in  [35],  and  a  third-order  BDF  scheme.  Concerning  accuracy,  we  have  first 
order  for  the  0-method,  second  order  for  the  leapfrog  and  BDF2  schemes,  and  third 
order  for  BDF3.  Concerning  stability,  the  amplification  factors  for  the  model  equation 

dq 

dt  =  ^ 

where  i  =  \f—l  and  w  e  R,  are  plotted  in  Figure  3.1  for  both  the  explicit  version 
(left)  and  the  implicit  version  (right)  of  the  considered  time  integrators.  The  meth¬ 
ods  which  guarantee  unconditional  stability  are  the  implicit  TM,  LF2,  and  BDF2. 
More  details  can  be  found  in  [33].  Having  completely  defined  our  time-integration 
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Fig.  3.1.  Amplification  factors  for  the  time  integration  schemes  listed  in  Table  3.1  for  the 
explicit  (left)  and  the  implicit  (right)  versions.  For  the  TM  and  LF2  methods  the  choice  6  =  0.5  has 
been  considered. 


method,  it  is  important  to  notice  that,  in  most  numerical  weather  prediction  codes,  the 
semi-implicit  method  is  usually  implemented  as  an  implicit  correction  to  an  explicit 
predictor  substep  [43].  In  our  case,  this  leads  to  the  following  algorithm: 

(i)  compute  the  explicit  predictor  step 

2  2 

cT  =  2  +  7 At  ^  /3m^(q"-m); 

m=0  m= 0 

(ii)  let  pm  =  crm  Pm ,  with  /?_ i  =  0  and  compute  q*  =  p_iqea:  +Em=o  P™q"_m; 

(iii)  let  38  =  ( ^  —  7  A t«£?°)  ,  where  is  the  identity,  and  compute  the  implicit 

corrector  step 

~  P- 1  (7A  t&X))  q  tt  =  -  (7At^°)  (  j]  pmq”“m 

\m= 0 

+  p- 1  (yAt^f0)  ; 

j2 

notice  that  qtt  is  an  approximation  of  (see  [33]); 

(iv)  update  the  solution  by  setting  qn+1  =  ^-(q u  ~  SLo 

4.  Spatial  discretization.  In  this  section,  we  address  the  spatial  discretization 
of  (2.5)  by  resorting  to  a  high-order,  nodal,  DG  formulation.  The  general  framework 
for  such  discretization  is  provided  by  [6,  5],  while  we  refer  to  [30,  34,  32]  for  the  aspects 
specifically  related  to  the  high-order  approximation. 

4.1.  Notation.  Let  be  an  open  bounded  domain  of  R2,  with  boundary  dfl 
and  outward  unit  normal  vector  n^o,  and  let  3h  denote  a  partition  of  f l  into  A fei 
nonoverlapping  curvilinear  quadrilateral  elements  K  which  are  images  of  the  reference 
element  K  =  [—1,  l]2  under  smooth,  bijective  maps  & k 

y/t  e  ,%  :  K  =  &K  (A)  . 

The  diameter  of  K  is  hx  and  we  let  h  =  max^e^y  hx-  The  set  of  the  edges  e  of  the 
triangulation  is  denoted  by  $h-  Let  8K  and  ne  dx  denote  the  boundary  of  K  e  3h  and 
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the  outward  unit  normal  vector  on  each  edge  e  e  dK ,  respectively,  and  d K  and  ru 
denote  the  boundary  and  the  outward  unit  normal  vector  for  the  reference  element 
K.  The  notation  x  =  with  x  =  (x,  z)  e  fl,  £  =  (£,  C)  e  K  will  be  used. 

We  also  associate  with  each  local  map  the  Jacobian  Jk  =  and  the  determinant 
.7 k |  ■  Although  not  essential,  we  will  assume  that  3dh  is  conforming,  that  is,  given 
Ki,  K2  e  we  either  have  that  Ki  n  K2  is  empty,  or  it  is  a  vertex  or  a  complete 
edge  e  e  <4.  Following  the  usual  notation  in  the  context  of  DG  formulation,  we  now 
define  averages  and  jumps  (see  [16,  3]).  Letting,  thus,  Xh  denote  a  generic  function 
piecewise  continuous  on  74  5  for  a  given  element  K ,  edge  e  c  dK  and  point  x  e  e  we 
define 


Xh{xint(K))=  Hm  x/i(y)j  Xh{xeXt(K))=  jim  Xfc(y). 

y  ->•  x  y  ->  x 

y  e  K  y  $K 

This  definition  can  be  extended  to  vector-valued  functions  by  applying  it  component¬ 
wise.  For  e  =  dK  n  dK' ,  x  6  e,  the  average  and  jump  for  a  scalar  function  Xh  and 
vector  function  r /,  piecewise  continuous  on  74  are  defined,  respectively,  as  follows: 

{x/l}(x)  =  i(x.(x“tW)+Xft(x“t^'))) 

(41)  \rh}  (X)  =  i(r;i(x“^))  +  r^x™4^'))) 

M(x)  =  Xh&ntlK))ne,8K  +Xh(xint{K,))ne,dK' 

[[[r/jjlKx)  =  rh(xint(Kl)  ®  neidK  +  rh(x.int(K''>)  ®nefSK>. 

Notice  that  (4.1)4  differs  from  the  usual  definition  of  the  jump  for  a  vector- valued 
function,  and  this  latter  can  be  recovered  as  [[17,]]  =  Tr  ( |J[r /,]]]).  Finally,  we  let 

maxjx^x)}  =max{xh(xintiK)),Xh(x-intiK'))}- 

4.2.  High-order  polynomial  space.  The  logically  square  structure  of  the  ref¬ 
erence  element  K  significantly  simplifies  the  construction  of  high-order  polynomial 
bases,  since  the  multidimensional  basis  can  be  obtained  as  a  tensor-product  of  the 
one-dimensional  basis.  For  an  integer  k  4  1,  letting  Pfc([o,  b])  denote  the  space  of 
polynomial  functions  of  degree  less  than  or  equal  to  k  on  [a,  6] ,  we  introduce  the 
following  two  bases  for  Pfc([— 1,1]):  {y>i(£)}i=o  the  Lagrangian  (i.e.,  nodal)  basis 
associated  with  an  arbitrary  set  of  nodes  e  [—1, 1],  while  {^i(£)}j= 0  is  the  Legendre 
(i.e.,  modal)  basis  (see  [11]).  We  also  define 

&ij(£)  =  <Pi(0<Pj(0  and  =  j (^iJ  xeA 

(0,  x$K 

For  simplicity,  a  cumulative  index  /,  ranging  from  1  to  A f  =  (k  +  1)2A fei,  is  also  bi- 
univocally  associated  with  the  indexes  (K,ij).  Although,  in  principle,  several  choices 
for  the  nodes  £,  are  possible,  a  convenient  one  is  represented  by  the  Legendre-Gauss- 
Lobatto  (LGL)  points,  defined  as  the  roots  of  the  polynomial  (1  —  £2)^fc(£)-  As  a 
matter  of  fact,  these  points  are  endowed  with  a  Gaussian  quadrature  rule  that  can 
be  exploited  to  improve  the  efficiency  of  the  resulting  scheme  (see  [11]).  The  finite 
element  space  is  now  defined  as  14  =  span{$i},  I  =  1, ...  ,7V.  Notice  that  functions 
in  14  are  in  general  discontinuous  across  edges  e  e  <4- 
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As  it  will  be  clear  from  the  forthcoming  sections,  the  DG  formulation  requires 
the  evaluation  of  both  two-dimensional  and  one-dimensional  integrals  on  K  and  dK, 
respectively.  To  this  end,  following  [30]  the  following  approximate  quadrature  rules 
will  be  adopted 


k 

(4.2)  f  x(x)  dx  =  [  x(x(£))|d*r(£)|d£  w  X  x(x(&,  <j))|  J*(&, 

•N'  JK  i,j=0 


and 


(4.3) 


da 


f  x(a)da=  2  f  x(x(ct))  \JK(a)  I  4T(»)naa 
llK  k,:kJb 

k 

*  2x(x(6,Co))|JAte,Co)|  |^Tte,Co)ngii? 

i= 0 

k 

+  XI  x(x(6,0)) \Ji<((k,Cj)\  JxT^k,Q) nSaA 

3=0 

k 

+  X*(x(&’Cfc))  \Ji<(£,i,(k)\  |j^T(^,Cfc)ng3>(, 


Wi 


2  =  0 
k 


+  X  Mo,( j))  \Jk(M)\  Kk-T(^0) Ci)ng4il 

3=0 


I< 


dK 


Wj 


Wi 


u3> 


where  x  is  a  generic  function  piecewise  continuous  on  ^54,  |]  •  ||  2  is  the  Euclidean  norm 
of  a  vector,  4  and  Q  are  Legendre-Gauss-Lobatto  points  previously  introduced,  and 
Wi  are  the  associated  weights,  defined  as 

1 


k(k  +  1)  \ipk(& 

In  the  following,  for  i,  j  =  0, . . . ,  k,  we  will  let 

WK/ij  Jk  (G  ■  Cj  )  UyiCy  , 


Jr  (^U  C?  )ne,<5iC 

UK,ij  =  < 

0 

V. 

difT(^i)  Cj)ne,8K 

(j,e)  e  {(0,  ei),  (fc,  e3)} 
(be)  6  {(0,e2),(fc,e4)} 
otherwise. 


Notice  that  the  degree  of  exactness  of  the  quadrature  rules  (4.2)  and  (4.3)  is  2 k  — 
1.  Nonetheless,  as  discussed  in  [30],  this  approximation  does  not  spoil  the  order  of 
accuracy  of  the  resulting  DG  formulation. 


4.3.  Discontinuous  Galerkin  discretization.  We  first  rewrite  (2.5)  introduc¬ 
ing  the  auxiliary  variable  S  as: 


(4.4) 


V- Te(q)  +  V  •  -Fu(q,  5)  +  G(q)  -  Tcv  (q  -  qb) 
5  —  Vq  =  0 


to  be  solved  in  f l  x  (0,Tfin]  with  suitable  initial  and  boundary  conditions.  An  ap¬ 
proximation  (qh,Sh)  =  (qb(x,  t),iSb(x,  t))  to  the  solution  (q(x,  t),  5(x,  t))  of  (4.4)  is 
sought,  such  that  (q h,Sh)  e  ( 14 )4  x  ( 14 )8  at  each  time  level.  In  the  following,  for 
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the  sake  of  clarity  we  will  often  omit  the  dependence  of  qb  and  Sh  on  (x,  t).  Also, 
all  products,  differential  operators,  and  average  and  jump  operators  are  intended  to 
be  applied  separately  to  each  density,  momentum,  and  energy  component  of  their 
arguments.  Multiplying  the  two  equations  (4.4)  by  test  functions  Xh  6  (Vk)4  and 
7 Zh  e  (I4)8,  respectively,  integrating  over  K  e  S'h  and  replacing  the  exact  solution  by 
its  approximation,  we  have 

^  J  q hXh  dx  =  -  J  V  •  [Te(<lh)  -  Fv{qh,  <Sft)]  Xh  dx 

+  f  [G(q;, )  -  Tcv  (qh  -  qb)]  Xh  dx,  MXh  e  (K)4 

Jk 

f  Sh  •  7Zh  dx  —  f  Vqh  •  7Zh  dx  =  0,  V7lhe(Vh)8. 

JK  JK 


Then,  formally  integrating  by  parts  and  introducing  the  numerical  fluxes  J~e  (qb ) , 

q(cth,Sh),  and  5(q;,,  Sh),  we  obtain  the  following  discrete  problem: 

find  (q/i(-,f),  Sh{-,t))  e  (14)4  x  (V/,)8  such  that  for  all  t  e  (0,  Tfin\  and  VAT  e 

(4.5) 

^  J  <lhXhdx  =  J  [^(q h)  -^(q h,Sh)]  ■  Vxbdx 

-f  [^e(q„)-^(q,5)l 

JdK  L 


n 8KXh  dcr 


lib 


+  [G(q h)  ~  Tcv  (qh  -  q6)]  Xhdx,  V xh  e  (Vh)4 

Jk 


l 


Sh  ■  TZh  dx  +  qb  V  •  7 Zh  dx  - 


qn,-/<  •  Kh  dcr  =  0, 

JdK 


V7 Zh  e  ( Vhf . 


Specification  of  the  numerical  fluxes  completes  now  the  definition  of  the  scheme. 
Concerning  the  hyperbolic  flux,  we  consider  here  the  Rusanov  flux  (see  [45,  54]) 


(4.6) 


where  |A|  =  max{ |v/(  -n|  +ab},  vb  and  ab  denoting  the  wind  velocity  and  sound  speed 
associated  with  qb,  respectively.  The  sound  speed  is  in  turn  defined  as  a  =  \JjRT, 
with  7  =  cp/cv.  We  notice  that  more  sophisticated  choices  than  (4.6)  are  also  possible, 
such  as  the  HLL,  HLLC,  or  Roe’s  fluxes  (see  [54]  and  also  [19]  for  comments  about 
the  effect  of  the  numerical  flux  on  the  computed  solution).  Concerning  the  viscous 
terms,  the  Bassi  and  Rebay  method  of  [5],  which  is  a  particular  case  of  the  local 
discontinuous  Galerkin  approach  described  in  [17],  is  adopted.  We,  thus,  set 


(4.7) 


q={q,,},  S={Sh}. 


Remark  4.1.  When  Tcv  vanishes,  i.e.,  far  from  the  open  boundaries,  the  con¬ 
servation  properties  of  (4.5)  follow  naturally  from  its  flux-form,  which  can  be  made 
explicit  by  taking  the  test  function  Xh  equal  to  the  characteristic  function  of  K  e 
yielding 


d_ 

dt 


I  [j 


Fe(qh)  ~  JFv(q,S)  •n5/fdcr+  G(qb)dx, 


L' 


and  by  noting  that  the  numerical  fluxes  are  single  valued  on  Sh-  Notice  also  that 
discrete  conservation  is  guaranteed  regardless  of  the  approximate  evaluation  of  the 
area  and  boundary  integrals. 
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We  now  take  the  test  function  Xh  *n  (4.5)i  equal  to  4>/.  The  left-hand  side  can 
be  written  as 

U8)  4£qAdx-M^, 

with 

(4.9)  Mij  =  4>/4>./  dx  =  widij. 

JK 

The  right-hand  side  defines  the  operator 


(4.10) 


(j^s(q;i))/=  f  \Te ( q/j )  —  Tv ( q/, , ) ] 

JK 


f  [^*<q' 

JPK  L 


A)-^(q,5) 


V4>/  dx 


n8K$i  du  + 


f  G(q*)$ 

JK 


i  dx 


for  I  =  1,...,A/;  is  the  discrete  counterpart  of  SfiNS  in  (2.5).  Finally,  the 

remaining  terms,  associated  with  the  damping  layer,  will  be  considered  in  section  4.4. 

4.4.  The  linear  problem.  Proceeding  as  done  in  section  2.2  to  define  ££NS 
in  (2.7)  from  S^NS  in  (2.5),  we  identify  within  (4.10)  the  following  linear  operator 


(4.11)  (£?hSqh)I  =  f  ^(q/i)-V4>/dx—  f  (qh)-ndK$i  da+ 
JK  JdK 

for  I  =  1, . . . ,  A f.  In  (4.11)  the  linear  numerical  flux  is 


f  G(q/l)$/clx, 
JK 


-  {r*}  +  ^  W, 

where  |A|  =  a  =  ■yJ'yRT  and  T  is  the  temperature  of  the  reference  state,  which  is 
assumed  to  be  a  continuous  function  over  O.  We  notice  that  the  linear  operator 
represents  the  approximate  counterpart,  using  the  DG  finite  element  method  supplied 
with  the  Rusanov  numerical  flux,  of  the  corresponding  space-continuous  operator 
2 '~£NS  introduced  in  (2.7).  In  other  words,  represents  in  itself  a  consistent  DG 

discretization  of  the  linearized  Euler  equations.  In  addition,  in  order  to  deal  with  the 
damping  terms  in  (4.5)i,  we  define 

(4.12)  (J^V),  =  -  f  T^q^/dx,  f°  =  \  T-q^/clx, 

JK  JK 

for  I  =  1, . . . ,  A f. 

5.  The  fully  discrete  problem.  The  fully  discrete  space-time  approximation  of 
problem  (2.5)  is  obtained  by  properly  substituting  into  the  multistep  time-advancing 
algorithm  illustrated  in  section  3  the  time  derivative  (4.8)  and  the  discrete  operators 
y^s{(\h),  Af);V'Sq/i,  _£?°q h  and  f°  introduced  in  (4.10),  (4.11),  and  (4.12),  respectively. 
At  each  time  level,  Sh  is  computed  from  (4.5)2- 

5.1.  Filtering  the  high-frequency  modes.  The  DG  method  using  high-order 
basis  functions  can  be  regarded  as  a  spectral  element  method  with  no  continuity 
constraint  among  neighboring  elements.  Since  high-orcler  methods  do  not  present  any 
intrinsic  numerical  diffusion,  they  are  prone  to  instabilities  due  to  nonlinear  mixing 
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and  the  Gibbs  phenomenon,  particularly  in  the  case  of  poorly  resolved  flows  (see  [9]). 
The  usual  way  of  dealing  with  this  instability  in  the  context  of  spectral  formulations 
is  to  introduce  a  filtering  operator  which  damps  the  high  frequency  modes  without 
altering  the  low  frequency  modes.  This  is  done  by  transforming  from  nodal  represen¬ 
tation  to  modal  representation,  applying  a  lowpass  filter,  and  then  by  transforming 
back  to  nodal  representation.  In  the  present  paper,  filtering  is  performed  using  a 
similar  strategy  as  in  [30],  and  the  action  of  the  filter  is  included  in  .  It  should 
be  mentioned  that  the  DG  need  for  a  filter  arises  only  due  to  our  particular  choice 
of  inexact  integration  rules  (of  order  2k  —  1);  choosing  2k  integration  rules  obviates 
the  need  for  a  filter,  but  the  cost  of  the  method  increases  due  to  the  need  to  invert 
a  nondiagonal  (albeit  small)  mass  matrix.  Furthermore,  having  a  nondiagonal  mass 
matrix  complicates  the  construction  of  the  pseudo-Helmholtz  problem;  this,  in  fact, 
is  the  key  to  the  success  of  the  semi-implicit  time-integration  method. 


5.2.  Solution  algorithm  for  the  pseudo-Helmholtz  operator  problem. 

The  linear  system  arising  from  the  SI  time  discretization  is  usually  dealt  with  by 
properly  combining  the  continuity,  momentum,  and  energy  equations,  in  such  a  way 
to  obtain  an  algebraically  equivalent  problem  of  diffusion-reaction  type,  referred  to  as 
the  pseudo-Helmholtz  operator  problem,  for  the  sole  pressure  variable.  This  equiv¬ 
alent  reformulation  has  the  computational  advantage  of  reducing  considerably  the 
dimension  of  the  problem;  moreover,  if  an  iterative  solver  is  adopted,  it  produces  a 
significant  acceleration  of  the  convergence  rate  and  simplifies  the  definition  of  the 
stopping  criterion.  Based  on  these  considerations,  an  extension  of  the  above  ap¬ 
proach  to  the  present  semi-implicit  DG  setting  is  highly  desirable,  albeit  being  far 
from  trivial.  Nevertheless,  such  an  extension  can  be  made  possible  by  conveniently 
exploiting  the  structure  of  the  approximate  quadrature  rules  (4.2)  and  (4.3).  By  do¬ 
ing  so,  one  obtains  a  formulation  that  can  be  regarded  as  a  LDG  discretization  of  a 
diffusion-reaction  problem  for  the  pressure  variable  where  the  auxiliary  flux  unknown 
is  statically  condensed  out  by  proper  use  of  mass  lumping.  In  this  sense,  the  resulting 
discrete  scheme  shares  some  similarities  with  hybridized  dual  mixed  methods  where 
static  condensation  is  the  crucial  approach  to  obtaining  a  linear  algebraic  problem 
for  the  sole  primal  variable  (see  [4]  for  an  introduction  to  static  condensation  at  the 
discrete  level,  and  [15,  14]  for  a  more  recent  development  on  this  subject).  In  the 
following,  the  two  terms  static  condensation  and  pseudo-Helmholtz  form  will  be  used 
interchangeably. 

For  simplicity,  we  assume  here  that  periodic  boundary  conditions  are  prescribed 
and  we  do  not  include  the  gravity  terms  into  the  implicit  part  of  the  problem.  We 
set  V  =  [U,  W~\T  and  denote  by  Fjf,  and  Fjf  the  rows  in  corresponding 
to  density,  momentum,  and  energy,  respectively.  Under  these  assumptions,  the  linear 
problem  arising  from  the  SI-DG  formulation  reads: 

find  (ptt,ytt,Ett,ptt)  e  14  x  (Vh)2  x  Vh  x  Vh  such  that  MK  e^,/  =  l,..  .,Af: 

(5.1) 


Ptt^i  dx  —  a  Vtt  •  V<I>/  clx  +  a 


pu^i dx  -  a 
jk  Ji 

t$i  dx  -  a  J 


I 

1 


:  t 

.  K 


Ett&i  dx  —  a  hVtt  *  V$j  dx  +  a 


dK 


JdK 


Ff 

•  n,-/cT/ 

da  = 

[ 

p*4>/  dx 

'k 

da  = 

[ 

V*<F,  dx 

- 

]K 

Ff 

•  iV/Uf, 

da  = 

[ 

E  *i>,  dx 

4- 

Ptt®i 

dx  =  — 

f  Ett$r 

Ik 

c 

"V 

Jk 
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with  a  =  p-i^fAt.  Equations  (5. 1)2, 3,4  then  immediately  provide  the  problem: 
find  ( ptt ,  Vtt)  e  Vh  x  (Vh)2  such  that  MIC  e  &h,  I  =  1, . . . , TV" 


Ptt^i  dx  -  a- 


R  r  v. 


L 


cv  JK 


h 


hVtt  •  V$j  dx  +  a 


■  j 


Vtt  •  3>j  dx  —  a  •  <&/  dx  +  a 


j  /iV’iifiA'fjdff 

JdK 

(pndK)  -&ida 

JdK 


L 

L 


p*$>i  dx 
V*  •  $/dx, 


where  3>/  =  [<!>/,  0]T  or  3>/  =  [0,4>/]T,  p*  =  -^E*,  and  the  numerical  fluxes  are 


(5.3)  V  =  {Vtt}  +  p  =  {pu}l+  ^DTVttl. 

2  Rhh  2 

Problem  (5.2)  can  be  fully  regarded  as  the  LDG  discretization  of  the  following  elliptic 
problem  for  the  sole  pressure  variable  (see  equations  (2.2)  and  (2.3)  in  [13]) 

(5.4)  -a2 V  •  (a2Vptt)  +  Ptt  =  -aV  •  (a2V*)  +  p* 

supplied  with  numerical  fluxes  (5.3),  except  for  the  fact  that  [[[Vtt]]]  is  used  in  (5.3)2 
instead  of  [[Vtt]]. 

Remark  5.1.  The  elliptic  numerical  fluxes  (5.3)  descend  from  the  hyperbolic 
numerical  flux  (4.6).  In  particular,  the  dissipative  term  ^-[[q/i]]  in  this  latter  flux 
gives  rise  to  the  two  stabilization  terms  in  (5.3). 

As  discussed  in  [13],  the  inclusion  of  the  jump  term  in  (5.3)2  makes  it  impossible 
to  compute  Vtt  element  by  element  in  terms  of  ptt  from  (5.2)2-  This  implies  that, 
starting  from  (5.2),  it  is  not  possible  to  obtain  in  an  efficient  way  a  discrete  counterpart 
of  (5.4)  involving  the  sole  pressure  variable.  However,  we  show  now  how,  by  taking 
advantage  of  the  approximate  quadrature  rule,  it  is  possible  to  compute  Vtt  node 
by  node  in  terms  of  ptt-  To  this  end,  we  need  to  work  out  the  matrix  formulation 
of  (5.2),  and  some  additional  notation  is  required.  We  denote  by  M  the  number  of 
quadrature  nodes  in  C7h,  and  assume  without  loss  of  generality  that  the  first  A4^  of 
such  nodes  are  located  on  <#),,.  For  the  Qth  quadrature  node,  we  denote  by  Iq(Q)  the 
degrees  of  freedom  collocated  at  the  quadrature  node  itself,  with  q  =  1 ,uq.  For 
Q  >  A4S  we  have  uq  =  1;  for  Q  ^  A4^,  on  the  one  hand,  we  have  uq  ^  2  thanks 
to  the  periodic  boundary  conditions  and,  on  the  other  hand,  the  regularity  of  C7h 
implies  an  upper  bound  for  uq.  Symmetrically,  for  a  degree  of  freedom  I  we  denote 
by  Q(I)  the  corresponding  quadrature  node.  Notice  that,  for  a  continuous  function 
X,  we  have  Xiqi(Q)  =  Xin(Q )  =  XQ-  For  a  given  pair  (Q,e),  with  Q  collocated  on  e, 
we  use  the  shorthand  notation  Iq  £  ,  I'q  e  to  indicate  two  degrees  of  freedom  such  that 
/Q,e  #  I'q  e,  Q(lQ,e)  =  Q{I'q  ancl  lQ,e  and  I'q  e  belong  to  a  couple  of  elements  IC, 
K'  such  that  8IC  n  8K'  =  e.  Notice  that  the  subscript  (Q,  e)  will  be  usually  omitted, 
since  it  is  clear  from  the  context.  Finally,  for  I  =  (K,ij)  we  denote  by  Si  the  set  of 
edges  e  such  that  e  £  dK  and  Q{I)  is  collocated  on  e. 

The  integrals  on  the  interior  of  K  can  be  easily  expressed  in  terms  of  the  TV"  x  J\f 
matrices  M,  Dx,  and  Dz  defined,  respectively,  in  (4.9)  and  by 

(5.5)  (Dx)jj  =  -  £  £^<t>./  clx,  (Dz)jj  =  -  £  clx, 

which,  due  to  the  discontinuous  nature  of  the  basis  functions,  are  block-diagonal. 
Concerning  the  boundary  integrals,  we  illustrate  the  treatment  of  (5.2)i,  as  analogous 
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considerations  apply  to  (5.2)2-  Dropping  the  subscript  tt  for  simplicity,  we  have 


I  —  h  V  •  XlpK^I  do- 

JdK  cv 

=  k’f  |  ^7^  ( Pi  —  Pi1 )  +  2  — —  [( Ui  +  Up )  nx  +  ( Wj  +  Wp )  nz]  1 

eeSi  V  v  ) 

=  2  +  H  coeI[(Ui+Ur)nx  +  (WI +Wp)nz]. 

Z  eE<S?i  Z  Cv  ee<?j 

Let  now  Ds ,  BNX,  and  BNZ  denote  the  Af  x  Af  matrices  such  that 


(£>*%)/  =  X  E  w!(a/-a/'), 

eeSi 


[BNXjZqh)I 


\  2  w/(9/  +qi')nx,z 


with  qh  denoting  either  ph,  Uh ,  or  ID/,, .  It  is  easy  to  verify  that,  up  to  a  permutation  of 
the  unknowns,  Ds ,  BNX ,  and  BiVz  have  a  block-diagonal  structure  with  nonzero 
blocks  of  dimension  tiq ,  Q  =  1, . . . ,  A4<?,  respectively.  In  particular,  for  the  case  of  a 
quadrature  node  belonging  to  one  sole  edge  e  (i.e.,  not  corner  point),  we  have  uq  =  2 
and  the  2x2  blocks 


Dq  =  - 
Q  2 


1  -1 

-1  1 


and 


BNxn  =  - 
XQ  2 


%xe,Iz(Q) 


xe,IpQ) 

xe,/2(0) 


tu 


Q> 


=  2 


t'2:e,J1(Q) 

^e,J2(Q) 


°ze,I1(Q) 

lze,I2(Q ) 


a;, 


Q’ 


with  =  wf0  e  =  wf^e. 

Summarizing,  the  matrix  counterpart  of  (5.2)  after  proper  use  of  numerical  inte¬ 
gration  reads 


(5.6) 


f  (M  +  aADs)p  +  a  [(Dx  +  BNX)  A  U  +  (Dz  +  BNZ)  A  W ]  =  Mp* 
I  (M  +  aADs)  U  +  a  (Dx  +  BNx)p  =  MU* 

(M  +  aADs)W  +a(Dz+BNz)p  =  MW*, 


where  A  and  A  are  AT  x  Af  diagonal  matrices  defined  as 


A ij  =  A q(j)  Sij,  Au 


—hQ(i)Sij. 

Uy 


The  key  element  for  the  reformulation  of  (5.6)  in  terms  of  the  sole  unknown  p  relies 
on  an  efficient  computation  of  the  inverse  of  the  matrix  A1DG  =  M  +  aADs.  We  have 


(Mdg)  1  =  (I  +aAM~1Ds)  1  M~l  :=  T,M~l. 

Since  E_1  =  (I  +  aA M~1DS),  it  turns  out  that  E  has  the  same  block-diagonal 
structure  as  Ds ,  so  that  its  computation  is  straightforward  (see  Figure  5.1). 
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Fig.  5.1.  Representation  of  the  block- diagonal  structure  of  the  matrix  E.  A  lxl  block  is 
associated  with  internal  quadrature  nodes  (A)  such  as  Q and  we  have  cr^1  =  1.  Quadrature  nodes 
belonging  to  one  sole  edge  (o),  such  as  Q 2,  give  rise  to  2x2  blocks.  Finally,  nodes  belonging  to  two 
or  more  edges  (□  and  o),  such  as  Q 3  and  Q 4,  give  rise  to  uq  x  nQ  blocks. 


Upon  defining  DGG  =  (Mdg)~1(Dx  +  BNX ),  DGG  =  { Mdg)~1{Dz  +  BN, ), 
p  =  E p*,  U  =  E/7*,  and  W  =  E W*,  system  (5.6)  can  be  written  as 


(5.7) 


p  +  a  \dxgA  U  +  D°gA  W 
U  +  aDGGp  =  U 


=  P 


W  +  aDGGp  =  W. 


Substituting  (5. 7)2, 3  into  (5.7)  1  we  obtain 

(5.8)  p  -  a2  (D°GAD°G  +  D°gAD°g)  p=p-a  (fjxG  A  U  +  DGG  A  w)  . 

Problem  (5.8)  is  the  discrete  counterpart  of  (5.4).  The  advantages  of  solving  (5.8) 
instead  of  (5.6)  will  be  numerically  demonstrated  in  section  6.1. 

6.  Numerical  results.  In  this  section,  the  numerical  validation  of  the  proposed 
scheme  is  carried  out.  To  simplify  the  presentation,  rather  than  giving  a  detailed  de¬ 
scription  of  each  test  case  setup,  we  provide  references  to  some  classical  works  in  the 
literature.  In  addition,  a  comprehensive  overview  of  all  the  test  cases  can  be  found 
in  [31].  For  ease  of  comparison  with  the  literature,  all  the  results  are  converted  to 
primitive  variables:  7 r  is  the  Exner  pressure,  u  and  w  are  the  horizontal  and  verti¬ 
cal  velocities,  respectively,  and  9  is  the  potential  temperature.  Deviations  from  the 
background  atmosphere,  which  is  characterized  by  a  uniform  Brunt-Vaisala  frequency 
N  =  |  jL  (see  [25]),  are  displayed.  Following  [48],  the  reference  state  q  (see  section  2) 
is  isothermal  in  all  the  test  cases,  with  T  being  the  highest  temperature  in  the  initial 
condition.  The  local  Courant  number  and  advective  Courant  number  are  defined  as 
C  =  anc^  CacLv  =  ’  respectively,  where  Hlgl  denotes  the  (variable)  spacing 

between  the  Legendre-Gauss-Lobatto  points  (see  section  4.2),  a  is  the  sound  speed, 
and  v  is  the  wind  velocity.  The  maximum  values  of  C  and  Cadv  in  the  domain  are 
denoted  by  (7max  and  G^ax,  respectively.  All  the  computations  described  in  this  sec¬ 
tion  employ  the  BDF2  scheme  of  Table  3.1;  notice,  however,  that  similar  results  are 
obtained  with  the  TM  and  LF2  schemes.  To  compare  our  results  with  other  results 
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in  the  literature,  it  is  convenient,  for  a  given  grid  of  high-order  elements,  to  define  the 
equivalent  resolution  as  the  resolution  of  a  uniform  grid  with  the  same  number  A4  of 
grid  points.  The  static  condensation  procedure  described  in  section  5.2  is  tested  for 
the  case  of  periodic  and  no-flux  boundary  conditions  and  implicit  treatment  of  the 
acoustic  waves  in  section  6.1.  For  the  remaining  test  cases,  the  linear  system,  which 
includes  an  implicit  treatment  of  the  gravity  waves,  is  solved  in  its  complete  form;  we 
plan  to  extend  the  static  condensation  procedure  to  the  general  case  in  future  work. 
Finally,  for  the  mountain  waves  test  cases  the  root-mean-square  (RMS)  errors  are 
computed  as  in  [31]. 

6.1.  Bubble  convection  experiments.  In  this  section,  we  study  four  idealized 
test  cases  characterized  by  buoyancy  driven  flows.  In  these  tests,  a  basic-state  atmo¬ 
sphere  is  considered,  which  is  assumed  to  be  at  rest  and  in  hydrostatic  equilibrium, 
and  a  thermal  anomaly,  with  a  consequent  density  perturbation,  is  introduced. 

The  first  test  case  aims  at  verifying  the  advantages  of  the  semi-implicit  time  dis¬ 
cretization  in  the  case  where  only  periodic  boundary  conditions  are  prescribed.  The 
computational  domain  is  the  rectangle  [Om,  1000  m]  x  [0  m,  2000  m],  and  the  initial 
datum  is  represented  by  a  thermal  anomaly  introduced  in  an  isothermal  atmosphere 
at  T  =  303  K  (notice  that,  thanks  to  this  choice,  the  deviations  from  the  reference  at¬ 
mosphere  are  zero  far  from  the  thermal  anomaly,  which  satisfies  the  periodic  boundary 
conditions  at  the  top  and  bottom  boundaries).  The  thermal  anomaly  has  amplitude 
—  15  A"  and  has  the  same  smooth  profile  described  in  [44].  Viscosity  is  set  equal 
to  0.2  m2/s.  The  computational  grid  is  composed  of  25  x  50,  10th-order  elements, 
with  equivalent  resolution  4  m,  while  the  time-step  is  0.02  s,  yielding  Cmax  =  6.1 
and  C™,ax  =  0.15.  The  computed  potential  temperature  is  displayed  in  Figure  6.1. 
For  this  particular  computation,  a  comparison  between  the  solutions  computed  with 
and  without  performing  the  static  condensation  shows  that  the  number  of  GMRES 
iterations  required  for  the  solution  of  the  linear  system  decreases  by  approximately 
a  factor  2,  while  a  rough  comparison  with  the  explicit  time  integration  indicates  a 
reduction  in  the  overall  computational  time  of  approximately  a  factor  5.  To  confirm 
this  result,  the  same  test  case  has  been  repeated  with  different  values  of  (7max,  using 
a  coarser  mesh  of  12  x  24,  10th-order  elements.  In  Figure  6.2  the  number  of  GMRES 
iterations  per  time  step,  the  total  number  of  GMRES  iterations,  and  total  CPU  time 


Fig.  6.1.  Cold  bubble  test  case,  potential  temperature  perturbation  at  time  levels  100  s  (loft) 
and  200  s  (right).  Contour  plots  using  values  between  —11.5  K  and  8  K  with  an  interval  of  1.625  K. 
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c  c 

Fig.  6.2.  Cold  bubble  test  case,  number  of  GMRES  iterations  per  time  step  (•,  left  axis),  total 
number  of  GMRES  iterations  (+,  first  right  axis),  and  total  CPU  time  (it,  second  right  axis)  as 
functions  o/ Cmax  for  the  case  of  static  condensation  (left)  and  no  static  condensation  (right). 


are  plotted  as  functions  of  Cmax  for  the  case  of  static  condensation  (left)  and  no  static 
condensation  (right) .  The  total  time  for  the  same  computation  using  an  explicit  BDF 
method  is  1.33  •  104  s;  using  the  largest  Courant  number  possible  we  note  in  Figure  6.2 
that  the  semi-implicit  method  with  static  condensation  yields  a  CPU  time  of  less  than 
2500  s,  which  is  a  factor  of  5  smaller  than  the  explicit  result.  From  these  results,  it  is 
clear  that  the  semi-implicit  time  discretization  increases  the  efficiency  of  the  scheme 
for  a  wide  range  of  Courant  numbers,  and  that  there  is  a  substantial  benefit  in  us¬ 
ing  the  static  condensation  procedure.  The  improvement  associated  with  the  static 
condensation  is  due  partly  to  a  reduction  in  the  number  of  GMRES  iterations  and 
partly  to  a  minor  cost  of  each  iteration,  since  the  unknown  is  represented  by  the  sole 
pressure  variable. 

The  second  test  case  is  similar  to  the  smooth  bubble  test  proposed  in  [12,  44]. 
The  basic  state  atmosphere  is  characterized  by  neutral  stratification,  and  the  flow  is 
driven  by  a  smooth  thermal  anomaly  of  which  the  maximum  amplitude  is  +0.5  K. 
Reflecting  boundary  conditions  are  applied  and  the  flow  is  inviscid.  The  compu¬ 
tational  domain  is  [0  m,  1000  m]  x  [0  m,  1000  m]  and  a  grid  composed  of  20  x  20, 
10th-order  elements  is  adopted,  with  equivalent  resolution  5  m.  The  time-step  is 
0.08  s,  yielding  Cmax  =  19  and  C^ax  =  0.12.  Figure  6.3  shows  the  computed  po¬ 
tential  temperature  perturbation  at  =  600  s  and  the  one-dimensional  profile 

along  z  =  700  m  for  the  same  quantity.  Minimum  and  maximum  values  for  the 
computed  solution  are  (—1.125  ■  10-5,  5.016  ■  10-6)  for  7r,  (—2.161  m/s,  2.161  m/s)  for 
u,  (-1.967  m/s,  2.758  m/s)  for  w,  and  (-7.303  •  10"2  if,  5.259  •  lO^/v)  for  0.  All 
these  results  are  in  good  agreement  with  those  reported  in  [31]  for  the  explicit  case. 
Concerning  conservation,  for  this  problem  we  expect  mass  and  total  energy  to  re¬ 
main  constant.  This  is  verified  up  to  machine  precision  in  the  numerical  simulation, 
where  we  observe  for  these  quantities  relative  deviations  equal  to  8.755  •  10-11  and 
4.627  •  10-11,  respectively. 

The  third  test  case  is  analogous  to  the  second  one,  except  that  a  nonsmooth  ther¬ 
mal  anomaly  is  considered  and  the  computational  domain  is  larger.  A  very  similar 
test  case  was  also  proposed  in  [44].  To  deal  with  the  nonsmooth  initial  datum,  a 
viscosity  v  =  0.4  m2 /s  is  introduced,  and  we  assume  Pi'  =  1.  Reflecting  boundary 
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Fig.  6.3.  Smooth  bubble  test  case,  potential  temperature  perturbation  at  time  level  600  s.  Left: 
contour  plot  using  values  between  0.025  A  and  0.525  .A  with  an  interval  of  0.05  K.  Right:  profile 
along  constant  height  z  =  700  m. 


x  x 

Fig.  6.4.  Nonsmooth  bubble  test  case,  potential  temperature  perturbation  at  time  level  600  s. 
Left:  contour  plot  using  values  between  0  K  and  0.52  K  with  an  interval  of  0.05  K.  Right:  profile 
along  constant  height  z  =  1000  m. 

conditions  are  applied,  except  for  the  energy  balance  equation,  where  the  temper¬ 
ature  gradient  is  imposed  on  bottom  and  top  boundaries  to  avoid  the  formation 
of  a  thermal  boundary  layer,  as  discussed  in  [31].  The  computational  domain  is 
[0  to,  1000  to]  x  [0  to,  1500  to]  and  a  grid  composed  of  20  x  30,  10th-order  elements  is 
adopted,  with  equivalent  resolution  5  to.  The  time-step  is  0.08  s,  yielding  C'max  =  19 
and  CZT  =  0.15.  Figure  6.4  shows  the  computed  potential  temperature  pertur¬ 
bation  after  600  s  and  the  one-dimensional  profile  along  z  =  1000  to  for  the  same 
quantity.  Minimum  and  maximum  values  for  the  computed  solution  are  (—2.563  • 
10-5, 1.418  •  10-5)  for  7 r,  (—2.386  m/s,  2.386  m/s)  for  u ,  (—3.965  7n/s,  2.450  to/s)  for 
w,  and  (—1.370  ■  10-3  K,  5.026  •  10_1  K )  for  9.  All  these  results  are  in  good  agreement 
with  those  obtained  for  the  same  problem  using  the  explicit  spectral  element  and  DG 
methods  described  in  [31].  Concerning  conservation,  for  this  problem  we  expect  mass 
to  remain  constant.  This  is  verified  up  to  machine  precision  in  the  numerical  simu¬ 
lation,  where  we  observe  a  relative  deviation  of  2.555  ■  10-11.  In  order  to  verify  that 
the  efficiency  gain  associated  with  the  semi-implicit  time  discretization  are  retained 
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Fig.  6.5.  Nonsmooth  bubble  test  case,  number  of  GMRES  iterations  per  time  step  (•,  left  axis), 
total  number  of  GMRES  iterations  (*,  first  right  axis),  and  total  CPU  time  (-k,  second  right  axis) 
as  functions  of  C'max  for  the  case  of  static  condensation  (left)  and  no  static  condensation  (right). 


3000  6000  9000  12000  15000  18000 


Fig.  6.6.  Density  current  test  case,  potential  temperature  perturbation  at  time  level  900  s  using 
400m  (left)  and  100m  (right)  resolution.  Contour  intervals  are  IK,  with  values  between  OK  and 
—  10 K ,  as  in  [51]. 

in  the  case  of  no-flux  boundary  conditions  and  nonsmooth  solutions,  the  number  of 
GMRES  iterations  per  time  step,  the  total  number  of  iterations,  and  the  total  CPU 
time  are  plotted  in  Figure  6.5  for  the  case  of  static  condensation  (left)  and  no  static 
condensation  (right),  for  10  x  15,  10th-order  elements.  The  total  time  for  the  same 
computation  using  an  explicit  BDF  method  is  2.12  •  104s;  using  the  largest  Courant 
number  we  note  in  Figure  6.5  that  the  semi-implicit  method  with  static  condensation 
yields  a  CPU  time  of  less  than  6000  s.  These  results,  analogous  to  those  of  Figure  6.2, 
show  the  robustness  of  the  semi-implicit  method. 

The  fourth  test  case  is  the  density  current  test  proposed  in  [51],  consisting  in  a 
cold  bubble  placed  in  a  neutral  atmosphere.  The  bubble  sinks  until  hitting  the  bottom 
boundary,  where  no-flux  conditions  are  imposed,  and  subsequently,  Kelvin-Helmholtz 
rotors  develop.  Viscosity  is  prescribed  in  such  a  way  that  a  grid-converged  solution 
can  be  obtained  at  approximately  50  m  resolution.  Figure  6.6  shows  the  computed 
solution  on  two  grids  composed  of  8  x  2  and  32  x  8,  8th-order  elements,  respectively, 
with  equivalent  resolution  400  m  and  100  in.  Notice  that,  thanks  to  the  symmetry 
of  the  problem,  the  solution  is  computed  only  in  half  of  the  domain.  The  time-steps 
are  0.8  s  and  0.2  s,  respectively,  with  Cmax  =  2.1  and  CU)ax  =  0.18.  The  first  case  is 
representative  of  a  poorly  resolved  flow,  since  the  resolution  is  too  coarse  to  capture  all 
the  features  of  the  grid-converged  solution,  while  the  second  case  is  representative  of  a 
well-resolved  flow.  By  comparing  the  results  in  Figure  6.6  with  those  in  [51],  it  can  be 
seen  that,  on  the  one  hand,  in  the  poorly  resolved  case  one  of  the  three  rolls  present  in 
the  reference  solution  is  clearly  reproduced,  yielding  a  result  which  is  comparable  to 
the  low  order  monotonic  upstream  method  (MUPL)  and  piecewise  parabolic  method 
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Fig.  6.7.  Self-convergence  experiment  for  the  density  current  test  as  in  [51],  Figure  4;  L 2 
errors  in  °C  for  9  as  functions  of  the  space  and  time  resolutions  (axes  in  m  and  s,  respectively). 
The  space  refinement  curve  (O)  is  obtained  at  fixed  time  step  At  =  1.5625  ■  10— 3  s,  with  a  reference 
solution  computed  at  25  m  space  resolution.  The  time  refinement  curves  are  computed  at  fixed  space 
resolution  200m  (*).  100m  (*),  and  50m  (•),  with  a  reference  solution  at  At  =  1.5625  ■  10— 3  s. 


(PPM)  solutions.  This  is  significant  because,  on  the  contrary,  the  high-order  fully 
local  spectral  (FLS)  method,  in  this  case,  produces  completely  meaningless  results. 
On  the  other  hand,  in  the  well-resolved  case  the  solution  obtained  with  the  SI-DG 
formulation  is  similar  to  the  one  obtained  with  the  FLS  and  close  to  the  grid-converged 
solution,  being,  thus,  superior  to  those  obtained  with  the  MUPL  and  PPM  schemes. 
Concerning  conservation,  we  observe  a  mass  relative  difference  of  4.600  x  10-12  and 
1.818  x  10“ 11  in  the  coarse  and  fine  resolution  cases,  respectively.  To  provide  a  more 
quantitative  measure  of  the  accuracy  of  the  SI-DG  formulation,  we  carry  out  a  self¬ 
convergence  study  analogous  to  the  one  described  in  section  3  of  [51].  For  this  test,  a 
reference  solution  at  spatial  resolution  25  m  and  At  =  1.5625T0-3  s  is  used  to  estimate 
the  separate  effects  on  the  error  of  space  and  time  refinements.  The  space  refinement 
is  obtained  by  fixing  the  time  step  as  in  the  reference  solution  and  considering  the 
four  spatial  resolutions  400  m,  200  m,  100  m,  and  50  m,  while  for  the  time  refinement 
the  spatial  resolution  is  fixed  and  decreasing  time  steps  are  considered,  starting  from 
the  largest  stable  time  step  allowed  by  the  semi-implicit  time  discretization.  The 
resulting  L 2  norm  of  the  error  in  d1  normalized  by  the  square  root  of  the  domain  area, 
is  plotted  in  Figure  6.7,  which  can  be  directly  compared  with  Figure  4  of  [51].  The 
choice  of  the  two  horizontal  axes  is  such  that  the  point  of  the  time  refinement  curve 
at  200  m  with  At  =  0.4  s,  which  is  the  maximum  stable  time  step  at  this  resolution, 
has  the  same  abscissa  as  the  point  at  200  m  resolution  in  the  space  refinement  curve, 
and  so  on  for  the  other  two  time  refinement  curves.  This  allows  for  an  immediate 
comparison  of  the  errors  resulting  from  the  time  and  space  discretizations,  which  can 
be  read  from  the  space  refinement  curve  and  the  corresponding  time  refinement  curve 
at  the  same  abscissa.  Concerning  the  space  refinement  curve,  we  observe  that,  at 
coarse  resolution,  the  error  of  the  SI-DG  method  is  smaller  than  the  error  reported 
in  [51]  by  a  factor  2,  a  ratio  which  reaches  a  factor  10  at  the  highest  resolution. 
The  experimental  convergence  rates  are  0.66,  2.90,  and  4.78,  which  are  consistent 
with  a  convergence  to  the  theoretical  value  9.  Concerning  the  time  refinement  curve, 
we  observe  the  theoretical  convergence  rate  2.  Finally,  concerning  the  ratio  between 
space  and  time  discretization  errors,  we  observe  that,  for  resolutions  coarser  or  equal  to 
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Fig.  6.8.  Inertia- gravity  wave  test  case,  potential  temperature  at  time  level  3000  s  (coordinates 
in  km).  Left:  contour  plot  with  contour  values  between  —1.5  ■  10— 3  K  and  3  ■  10— 3  K  with  contour 
interval  5  ■  10— 4  K  (negative  values:  dashed  lines).  Right:  profile  along  5  km  height  for  the  semi- 
implicit  (continuous)  and  for  the  explicit  (dashed)  DG  formulations. 


100  to,  the  time  discretization  error  at  the  largest  stable  time  step  is  roughly  one  order 
of  magnitude  smaller  than  the  space  discretization  error,  so  that  the  larger  stability 
of  the  semi-implicit  time  integrator  can  be  fully  exploited  without  compromising  the 
quality  of  the  results.  The  situation  is  different  for  the  well-resolved  case  at  50  to 
resolution,  where  the  time  discretization  error  is  larger  than  the  space  discretization 
one.  For  this  case,  one  should  either  use  a  smaller  time  step  or  resort  to  higher  order 
methods,  such  as  BDF3  in  Table  3.1. 

6.2.  Inertia-gravity  waves.  In  this  section,  we  study  the  nonhydrostatic 
inertia-gravity  wave  test  proposed  in  [49].  A  background  atmosphere  is  considered 
with  constant  Brunt-Vaisaila  frequency  N  =  0.01  s-1  and  uniform  horizontal  flow 
u  =  20m/ s.  The  computational  domain  is  [0  km,  300  km]  x  [0  fcm,  10  fern.],  with  no- 
flux  boundary  conditions  on  bottom  and  top  boundaries  and  periodic  conditions  on 
lateral  boundaries.  The  flow  is  inviscid.  The  initial  condition  is  represented  by  a 
thermal  anomaly  centered  at  (x,z)  =  (100  km,  5  km),  and  the  flow  is  simulated  un¬ 
til  Tfin  =  3000  s.  A  grid  composed  of  120  x  4,  10th-order  elements  is  employed, 
with  equivalent  resolution  0.25km.  The  time-step  is  Is,  yielding  Cmax  =  5.1  and 
CT  =  0.24.  The  computed  potential  temperature  9  is  shown  in  Figure  6.8.  These 
results  are  in  good  agreement  with  those  presented  in  [49,  1,  31].  To  assess  in  greater 
detail  the  effect  of  the  semi-implicit  time  discretization,  in  Figure  6.8,  right,  the 
computed  profile  from  [31]  with  explicit  time  stepping  is  also  reported.  It  can  be 
seen  that  the  two  curves  are  nearly  coincident.  Extrema  for  the  computed  solution 
are  (  —  7.128  •  10-7, 9.106  •  10-7)  for  7 r,  (  —  1.061  ■  10-2  m/s,  1.064  ■  10_2to/s)  for  u, 
(-2.402  ■  10“3to./s,  2.877  •  10“3to/s)  for  w,  and  (-1.511  ■  10“3  A",  2.806  ■  10-3  K)  for 
9.  The  corresponding  relative  differences  with  respect  to  the  explicit  case  are  46.9%, 
0.56%,  13.4%,  and  0.68%,  respectively,  thus,  confirming  that  the  distortion  caused  by 
the  semi-implicit  time  discretization  for  acoustic  modes,  which  is  evident  in  the  large 
pressure  difference,  does  not  affect  the  slow  modes.  Finally,  concerning  conservation, 
we  notice  that  for  this  problem  we  expect  mass,  horizontal  momentum,  and  total 
energy  to  remain  constant.  This  is  verified  up  to  machine  precision  in  the  numer¬ 
ical  simulation,  where  we  observe  for  these  three  quantities  relative  errors  equal  to 
1.669  •  10-8,  2.645  ■  10-7,  and  1.640  •  10-8,  respectively. 

6.3.  Mountain  wave  simulations.  In  this  section,  we  study  three  test  cases 
based  on  the  simulation  of  hydrostatic  and  nonhydrostatic  mountain  waves.  Besides 
assessing  the  overall  accuracy  of  the  proposed  numerical  scheme,  we  intend  with  these 
tests  to  verify  the  robustness  of  the  approach  discussed  in  section  2.1  and  section  3  to 
deal  with  open-boundary  conditions  in  the  framework  of  a  conservative  formulation 
for  the  flow  equations  and  semi-implicit  time  integration.  All  the  test  cases  consider 
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a  uniform  horizontal  flow  impinging  over  an  isolated  mountain.  The  background 
atmosphere  is  characterized  by  constant  Brunt-Vaisala  frequency;  no  flux  conditions 
are  imposed  on  the  bottom  boundary  while  open-boundary  conditions  are  imposed 
on  lateral  and  top  boundaries.  The  flow  is  inviscid  in  all  the  test  cases.  A  classical 
review  on  mountain  waves  is  [50],  while  we  refer  to  [36,  24,  40,  7,  46,  29,  28]  for  a 
detailed  description  of  the  test  cases  and  reference  solutions. 

For  the  first  test  case,  we  consider  an  isothermal  atmosphere  at  T  =  250K 
and  Brunt-Vaisaila  frequency  N  =  1.95  •  10_2s_1  and  a  uniform  horizontal  flow 
at  u  =  20  m/s.  The  computational  domain  is  [0  km,  240  km]  x  [0  km,  30  km],  and  the 
mountain  profile  is  defined  by  the  versiera  di  Agnesi 

(6-1)  hm{x)  = - 77 , 

*  +  (*?*) 

with  hmo  =  1  to,  xc  =  120  km.,  and  ac  =  10  km..  As  shown  in  [50],  this  choice 
of  parameters  results  in  hydrostatic  flow.  A  grid  composed  of  20  X  12,  10th-order 
elements  is  adopted,  with  equivalent  resolution  1.2  km.  in  the  horizontal  and  0.25  km. 
in  the  vertical.  The  time-step  is  3.5  s,  yielding  Cmax  =  18.6  and  G™)//  =  0.18.  The 
computed  solution  at  time  level  =  lOhrs  is  presented  in  Figure  6.9.  A  pseudo- 
analytic  solution  computed  with  Fourier  transform  techniques  is  also  represented.  In 
addition,  Figure  6.11,  left,  shows  the  computed  vertical  flux  of  horizontal  momentum 
(whose  magnitude  equals  the  drag  exerted  by  the  flow  on  the  obstacle),  at  various 
time  levels,  normalized  by  the  analytic  value 

(6.2)  Mh  =  osuNh2mo , 

with  ps  =  p(0)  denoting  the  surface  density.  In  general,  an  overall  agreement  of  the 
computed  solution  with  the  analytic  solution  can  be  observed.  The  RMS  errors  are 
1.44  x  10-7  for  7 r,  2.61  x  10-3  m/s  for  u,  7.44  x  10-5  m/s  for  w,  and  1.79  x  10-3  K 
for  0.  Concerning  the  normalized  momentum  flux,  its  value  is  close  to  1  far  from 
the  upper  boundary,  and  goes  to  zero  within  the  damping  layer.  This  fact  is  in  good 
agreement  with  the  analytic  results  and  also  confirms  that  a  steady  state  configuration, 
characterized  by  a  uniform  momentum  flux,  is  attained. 

For  the  second  test  case,  we  consider  a  constant  stability  atmosphere  with  Brunt- 
Vaisaila  frequency  N  =  1  ■  10-2  s_1  and  surface  temperature  Ts  =  280  A'  and  a  uni¬ 
form  horizontal  flow  with  u  =  10  m/s.  The  computational  domain  is  [0  km.,  144  km]  x 
[0  km,  30  km],  and  the  mountain  profile  is  given  by  (6.1)  with  hmo  =  1  to,  xc  =  72  km 
and  ac  =  1  km..  This  choice  of  parameters  results  in  nonhydrostatic  flow.  A  grid  com¬ 
posed  of  40  x  10,  10th-order  elements  is  adopted,  with  equivalent  resolution  0.36  km.  in 
the  horizontal  and  0.3  km  in  the  vertical.  The  time-step  is  2  s,  yielding  Gmax  =  8.39 
and  =  0.17.  The  computed  solution  at  time  level  Tfin  =  5hrs  is  presented  in 

Figure  6.10,  together  with  a  pseudo-analytic  solution.  Figure  6.11,  right,  shows  the 
computed  vertical  flux  of  horizontal  momentum  at  various  intermediate  time  levels, 
normalized  by  the  analytic  value 

Mnh  =  -0A57^psuNh2mo. 

In  general,  an  overall  agreement  of  the  computed  solution  with  the  analytic  solution 
can  be  observed,  although  some  incorrect  maxima  can  be  seen  downstream  of  the 
obstacle,  which  can  be  possibly  ascribed  to  unphysical  reflections  in  the  absorbing 
layer.  The  RMS  errors  are  2.28  x  10-8  for  n,  6.98  x  10-4  m/s  for  u,  2.59  x  10-4  m/s 
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Fig.  6.9.  Linear  hydrostatic  mountain  test  case,  computed  solution  at  time  level  10  hrs  (black, 
dashed  line  for  negative  values)  and  analytic  steady  state  solution  (gray),  coordinates  in  km.  Left: 
horizontal  velocity  u  with  contour  lines  between  —2.5  •  10— 2  m/s  and  2.5  •  10— 2  m/s  with  contour 


interval  5-10  05  m/s.  Right:  vertical  velocity  w 
5  •  10~3  m/s  with  contour  interval  5  •  10— 4  m/s. 


with  contour  lines  between  —5  •  10  a  m/s  and 
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Fig.  6.10.  Linear  nonhydro  static  mountain  test  case,  computed  solution  at  time  level  5  hrs 
(black,  dashed  line  for  negative  values)  and  analytic  steady  state  solution  (gray),  coordinates  in 
km.  Left:  horizontal  velocity  u  with  contour  lines  between  —2.5  •  10 ~2  m/s  and  2.5  •  10 ~2  m/s  with 
contour  interval  2.5  •  10~3  m/s.  Right:  vertical  velocity  w  with  contour  lines  between  —5  •  10“ 3  m/s 
and  5  •  10~3  m/s  with  contour  interval  5  •  10— 4  m/s. 


Fig.  6.11.  Mountain  wave  test  cases,  normalized  vertical  flux  of  horizontal  momentum.  Left: 
hydrostatic  case,  time  levels  2  hrs,  4  hrs,  6  hrs,  8  hrs,  and  10  hrs.  Right:  nonhydro  static  case,  time 
levels  1  hrs,  2  hrs,  3  hrs,  4  hrs,  and  5  hrs. 


for  w,  and  2.50  x  10 -4iC  for  0.  Concerning  the  normalized  momentum  flux,  the  same 
considerations  as  for  the  previous  case  apply. 

The  third  mountain  test  case  has  been  proposed  in  [46] ,  and  various  solutions  are 
available  in  the  literature,  as  for  instance,  [29,  28].  A  uniform  stability  background 
atmosphere  is  considered  with  Brunt- Vaisaila  frequency  N  =  1  •  10~2s_1,  surface 
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Fig.  6.12.  Schar  mountain  test  case,  computed  solution  at  time  level  10  hrs  (black,  dashed  line 
for  negative  values)  and  linearized  analytic  steady  state  solution  (gray),  coordinates  in  km.  Left:  hor¬ 
izontal  velocity  u  with  contour  lines  between  —2  m/s  and  2m/ s  with  contour  interval  0.2  m/s.  Right: 
vertical  velocity  w  with  contour  lines  between  —2  m/s  and  2m/ s  with  contour  interval  0.05  m/s. 


Fig.  6.13.  Schar  mountain  test  cases.  Left:  computational  grid,  element  boundaries  are  dis¬ 
played  in  black  and  LGL  points  are  located  at  the  intersections  of  the  gray  lines.  Right:  normalized 
vertical  flux  of  horizontal  momentum  at  time  levels  2  hrs,  4  hrs,  6  hrs,  8  hrs,  and  10  hrs. 

temperature  Ts  =  280  K,  and  a  uniform  horizontal  flow  with  u  =  10  m/s.  The 
computational  domain  is  [—25  km,  25  km\  x  [0  km,  21  km\,  and  the  mountain  profile  is 


with  hmo  =  250  m,  ac  =  5  km,  and  Ac  =  4  km.  As  pointed  out  in  [46],  the  orogra¬ 
phy  profile  forces  two  distinct  types  of  internal  waves:  a  large-scale  hydrostatic  wave 
characterized  by  deep  vertical  propagation  and  a  smaller-scale,  nonhydrostatic  wave 
characterized  by  rapid  decay  with  height.  A  grid  composed  of  20  x  10,  10th-order 
elements  is  adopted,  with  equivalent  resolution  0.25  km  in  the  horizontal  and  0.21  km 
in  the  vertical.  The  time-step  is  1.4  s,  yielding  Cmax  =  8.40  and  C^ax  =  0.20.  The 
computed  solution  at  time  level  Tfin  =  10  hrs  is  presented  in  Figure  6.12,  together 
with  a  pseudo-analytic  linear  solution.  Notice,  however,  that  due  to  the  nonnegligible 
height  of  the  mountain,  the  linear  solution  should  be  taken  as  a  qualitative  reference 
rather  than  the  “truth”.  A  good  agreement  of  the  numerical  solution  is  observed 
both  with  the  analytic  linear  solution  and  with  results  in  the  literature.  In  particular, 
notice  that  the  numerical  solution  correctly  reproduces  the  different  vertical  struc¬ 
tures  of  the  two  superimposed  waves,  with  the  small-scale  perturbation  exhibiting  the 
correct  decay  with  height.  The  RMS  errors  are  7.93  x  10-6  for  7r,  1.87  x  10-1  m/s  for 
u,  3.86  x  10_2to/s  for  w,  and  4.41  x  10-2  K  for  0.  Figure  6.13  shows  a  detail  of  the 
computational  grid,  left,  and  the  computed  vertical  flux  of  horizontal  momentum  at 
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various  intermediate  time  levels,  normalized  by  the  analytic  value  for  the  hydrostatic 
case  (6.2),  right.  The  momentum  flux  is  constant  along  the  vertical,  indicating  that 
a  steady  state  configuration  is  reached,  and  a  value  different  from  one  is  explained  by 
the  different  shape  of  the  mountain  with  respect  to  the  previous  test. 

7.  Conclusions.  In  the  present  article  we  have  proposed  a  formulation  in  which 
the  SI  time  integration  strategy  is  adopted  in  the  context  of  a  high-order  DG  spa¬ 
tial  discretization  for  the  solution  of  the  nonhydrostatic,  compressible  Navier-Stokes 
equation  for  atmospheric  flows.  The  main  reason  for  investigating  this  combination 
is  the  increase  of  the  efficiency  of  the  DG  method  when  applied  to  mesoscale  flows, 
and,  more  in  general,  to  low  Mach  number  compressible  flows.  A  critical  step  in  a  SI 
formulation  is  represented  by  the  solution  of  the  linear  system  for  the  implicit  part  of 
the  scheme.  We  have  shown  that  it  is  possible  to  reformulate  such  a  problem  in  terms 
of  a  pseudo-Helmholtz  operator,  and  that  the  resulting  discretization  fits  into  the 
LDG  framework  for  elliptic  problems.  We  have  also  described  how  to  sidestep  some 
well-known  difficulties  in  dealing  with  penalization  terms  in  the  numerical  fluxes  by 
exploiting  the  LGL  numerical  quadrature.  The  potential  benefits  of  this  approach 
have  then  been  demonstrated  with  some  classical  numerical  tests.  In  the  future,  we 
plan  to  implement  the  pseudo-Helmholtz  form  of  the  linear  system  for  the  case  of 
general  boundary  conditions,  develop  high-orcler  semi-implicit  methods  with  adaptive 
time-stepping,  explore  alternative  high-order  non-rcflective  boundary  conditions,  and 
extend  the  model  to  three  dimensions  to  investigate  the  effects  of  rotation.  A  fur¬ 
ther  possible  extension  is  represented  by  the  introduction  of  the  semi-Lagrangian  DG 
method  proposed  in  [41]  to  deal  with  the  stability  limit  associated  with  advection. 
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